TEMMETA is a pure python package that aims to facilitate basic processing of TEM data in a completely free way. The goal is to provide access to some of the most common operations with a single command, and to provide convenient visualisation widgets that work in jupyter notebooks. In the future I hope to improve the integration with hyperspy and make TEMMETA essentially a propper extension of hyperspy.
We set the plotting "back-end", this means where the plots will appear. With the back-end "notebook" the plots will appear embedded in this notebook. You could also use qt which will make pop-up windows.
%matplotlib notebook
Here we import the most important TEMMETA modules. With data_io
(imported here as dio) we can open and extract data objects from EMD files. It's possible that as the project evolves, the names of modules and packages will be changed. With image_filters
we get access to a few functions for altering images such as gaussian and median filters.
from temmeta import data_io as dio
from temmeta import image_filters as imf
from importlib import reload
reload(dio)
reload(imf)
<module 'temmeta.image_filters' from '/Users/nielscautaerts/Documents/PythonProjects/TEMMETA/temmeta/image_filters.py'>
We also import a number of other utilities for working more directly with the data
# numpy is for working with data arrays
import numpy as np
# pyplot is for making and manipulating plots
import matplotlib.pyplot as plt
# variable necessary to properly deal with paths
import os
here = os.path.abspath('')
print(here)
/Users/nielscautaerts/Documents/PythonProjects/TEMMETA/examples
The data with which the rest of this notebook works is not bundled together and must be downloaded. It is available through public link on owncloud and will be downloaded here using the wget utility.
save_folder = os.path.join(here, "data")
# Run this cell if you need to download the data
import wget
url1 = "https://owncloud.gwdg.de/index.php/s/kZ5BkdNJmcbJoFM/download" #17_ 20191205
url2 = "https://owncloud.gwdg.de/index.php/s/emUX6ksqDM6wTR1/download" #SuperX-HAADF
if not os.path.isdir(save_folder):
os.mkdir(save_folder)
wget.download(url1, out=os.path.join(save_folder, "hrtem.emd"))
wget.download(url2, out=os.path.join(save_folder, "stem-edx.emd"))
emd1_path = os.path.join(save_folder, "stem-edx.emd")
emd3_path = os.path.join(save_folder, "hrtem.emd")
Here we import the first emd file and open it as an EMDFile object. It is a h5py File object with some extensions. This will allow us easy access to the datasets.
emd1 = dio.EMDFile(emd1_path)
You can visualize the datasets inside the emd file as in the next cell. This reveals the datasets that can in principle be read in and converted to TEMMETA dataset objects for further processing.
Note: sometimes you will see additional "image" datasets that are actually derived from the spectrum data. These datasets will not be importable. Please import the spectrumstream and then create your own images instead (see later)
emd1.print_simple_structure()
------------------------------ Image datasets ------------------------------ Dataset number: 0, UUID: 6fdbde41eecc4375b45cd86bd2be17c0, Shape: 256x256, Frames: 240, Data type: uint16, Min:17381, Max:57476 --------------------------------------- SpectrumStream datasets --------------------------------------- Dataset number: 0, UUID: f5a4ba0965a5444b8c46cc420cf7fef0, Length: 15933353, Data type: uint16,
To extract the data from the dataset and create a nice easy object to work with, you need to pass in the dataset type and UUID, which uniquely identifies the dataset. You can get this from the table above. Here we extract the first dataset, representing an image stack. TEMMETA will try its best to recognize what type of dataset it is and create the appropriate one.
stack = emd1.get_dataset("Image", "6fdbde41eecc4375b45cd86bd2be17c0")
You can see that "stack" variable now contains a GeneralImageStack object.
type(stack)
temmeta.data_io.GeneralImageStack
We can create an interactive visualisation of the image stack with the following cell. This allows you to add a scalebar, check the different frames in the stack, and adjust contrast and brightness. I wanted to implement a "save" button also which opens up a dialog to save that frame, unfortunately this is not working at the moment, you will have to create and save frames manually
stack.plot_interactive()
We can select a frame of the stack or we can average all the frames with get_frame
and average
respectively
frame123 = stack.get_frame(123)
av = stack.average()
These are now individual images, objects that have a different set of tools than the stack
print(type(frame123))
print(type(av))
<class 'temmeta.data_io.GeneralImage'> <class 'temmeta.data_io.GeneralImage'>
One of the things you could do with an image dataset is make an intensity profile along a certain line defined by two points. You could first plot the image in an interactive way using plot_interactive()
and then hover the mouse over where you would like to have the points of your line. The coordinates of the mouse will appear in the lower-right corner. Write these down and use them as coordinates.
ax, im = av.plot(dpi=50)
# here we define our lines and also add it to the plot as an arrow for clarity
x1, y1 = (35, 101)
x2, y2 = (205, 101)
ax.arrow(x1, y1, x2-x1, y2-y1, color = (0., 1., 0.), head_width=5)
x3, y3 = (204, 200)
x4, y4 = (110, 20)
ax.arrow(x3, y3, x4-x3, y4-y3, color = (1., 0., 1.), head_width=5)
<matplotlib.patches.FancyArrow at 0x14c20adf0>
To create the line profile do the following. You could optionally set the w (width of the line in pixels) parameter. Pixels perpendicular to the line are added together. By default w=1. Here we create 4 profiles, also from the one frame so you can see the difference.
# the green arrow in the averaged image
profil1 = av.intensity_profile(x1, y1, x2, y2)
# the green arrow in the averaged image but with a higher width
profil1_b = av.intensity_profile(x1, y1, x2, y2, w=3, averaged=False)
# the green arrow in the one frame, will be much more noisy
profil1_c = frame123.intensity_profile(x1, y1, x2, y2)
# the magenta arrow in the averaged image
profil2 = av.intensity_profile(x3, y3, x4, y4)
We can plot the profiles with their "plot" function. This returns an axis object and a list of plot objects (only one element in this case). You can use the axis to pass it to the plot function of another profile to plot them on the same figure.
# plot returns the axis object
ax, prof1 = profil1.plot()
_, prof1b = profil1_b.plot(axis=ax)
_, prof1c = profil1_c.plot(axis=ax)
_, prof2 = profil2.plot(axis=ax)
# we can modify the properties of the plots, like give each a label as such
prof1[0].set_label("Profile 1")
prof1[0].set_color("Green")
prof1b[0].set_label("Profile 1b")
prof1b[0].set_color("Green")
prof1b[0].set_linestyle("-.")
prof1c[0].set_label("Profile 1c")
prof1c[0].set_color("Green")
prof1c[0].set_linestyle("--")
prof2[0].set_label("Profile 2")
prof2[0].set_color("Magenta")
# then we can autogenerate a legend
ax.legend()
# some layout optimization
ax.figure.tight_layout()
If we are not happy with the units, we can change them using set_scale
. For example, we might prefer units of nm in the plot instead of m. So we will use the current pixelsize
of the object and multiply by 1e9. We also set the unit to nm. This will change and overwrite the values in the metadata.
profil1.set_scale(profil1.pixelsize*1e9, "nm")
Replotting shows that it was effective.
profil1.plot()
(<AxesSubplot:xlabel='x (nm)', ylabel='Intensity (a.u.)'>, [<matplotlib.lines.Line2D at 0x13e087700>])
Again a note on metadata. You can check the metadata of all derivative datasets again with the metadata
attribute. Any derivative dataset stores the metadata of it's parent inside the parent_meta
key. The process to arrive at the child dataset is stored in process
. In this way the complete processing chain is kept inside the dataset metadata, as long as you perform operations on the object itself and not directly on the data. For the profile, we have a chain of two parents. The original_metadata
in the measured dataset is no longer stored in derivative datasets.
With this philosophy of traceability, any time an operation is performed on a dataset, it will create and return a new object. The data from parent objects is not overwritten. If you want to keep using the same variable, conserve memory, and perform operations on an object in place, most functions have the inplace
argument which is False
by default. If you set this to True
, data
and metadata
will be updated, but the metadata will also show the operation-chain.
print(profil1.metadata.process)
print(profil1.metadata.parent_meta.process)
Intensity line profile from (35, 101) to (205, 101) with width 1 Averaged all frames in stack
profil1.history
* Intensity line profile from (35, 101) to (205, 101) with width 1 | * Averaged all frames in stack | * Extracted dataset Image/6fdbde41eecc4375b45cd86bd2be17c0 from EMD file stem-edx.emd
Finally the intensity profile can be exported to an excel or a hyperspy dataset. Hyperspy takes a long time to load into memory so it is only loaded into memory once an export needs to take place. A filename can be optionally provided to save the object to a hyperspy file. A hyperspy file is also a hdf5 file, with a much more logical structure than Velox EMD files. You can import the file again later and continue working on it with hyperspy. All metadata related to the acquisition however is not stored in the hyperspy dataset.
profil1.to_excel("profile1.xlsx")
profil1_hs = profil1.to_hspy("profile1.hspy")
Importing hyperspy, this can take a while... /Users/nielscautaerts/opt/anaconda3/envs/devel/lib/python3.8/site-packages/pyUSID/viz/__init__.py:16: FutureWarning: Please use sidpy.viz.plot_utils instead of pyUSID.viz.plot_utils. pyUSID.plot_utils will be removed in a future release of pyUSID warn('Please use sidpy.viz.plot_utils instead of pyUSID.viz.plot_utils. ' WARNING:hyperspy_gui_traitsui:The nbAgg matplotlib backend is not compatible with the traitsui GUI elements. For more information, read http://hyperspy.readthedocs.io/en/stable/user_guide/getting_started.html#possible-warnings-when-importing-hyperspy. WARNING:hyperspy_gui_traitsui:The traitsui GUI elements are not available.
Overwrite 'profile1.hspy' (y/n)? y
type(profil1_hs)
hyperspy._signals.signal1d.Signal1D
Images are of course the bread and butter of microscopy. Here I try to show some of the things you can do with the GeneralImage
class. Intensity line profiles were already demonstrated before. We continue working on our averaged image calculated from the stack.
The main thing you will want to do with images is visualize them, apply filters, measure distances and angles, calculate Fourrier transforms, etc. For a quick look, we can also plot images with some interactive sliders. The sliders can at times be buggy which is an issue with matplotlib widgets. I recommend clicking inside the sliders instead of dragging them.
av.plot_interactive()
First we will also change the scale from m to nm on the image for convenience. While the scalebar is smart, for all kinds of measurements nm units will be more convenient. In addition, derived datasets, such as line profiles will also use the same units.
Note: call this cell only once. If you call it twice then it will look again what the current pixelsize is and again multiply by 1e9.
av.set_scale(av.pixelsize*1e9, "nm")
You can verify that it worked by checking the dimensions. this gives the size of the image in real units in a 4-tuple
av.dimensions
(4.508282741280623, 'nm', 4.508282741280623, 'nm')
Visualisation happens mainly through the plot command, it has innumerable options but they are a bit hidden. The main direct keyword arguments are width
which determines the plot size on the screen (the scalebar is not scaled with the image size), scalebar
which you set to True
(default) or False
to add a scalebar, show_fig
(T/F) to show the figure or not (if it's only for export), filename
(optional) for when you directly want to write out the plot to a file. You can pass arguments to the scalebar via sb_settings
which accepts a dictionary of settings, while imshow_kwargs
accepts a dictionary passed to the imshow
command in pyplot. If you want to modify how the scalebar looks sb_settings
is very important. Since the Axisimage object is returned, you can still change the settings of the image (like the colormap, the intensity scaling, etc) after you call the plot command.
For all the possible sb_settings
, check out the documentation of matplotlib-scalebar
For all the possible imshow_kwargs
, check out the documentation on imshow
Here I just show one example of a non-standard plot
scaleset = {"location": 'lower right',
"color": "white",
"length_fraction": 0.25,
"frameon": True,
"box_color": "black",
"font_properties": {"size": 14, "weight": "bold"}}
imgset = {"cmap": "jet", "vmin": 2e4, "vmax": 5e4}
ax, im = av.plot(width = 10, sb_settings = scaleset, imshow_kwargs = imgset)
You could always change the imshow
properties by accessing the Axesimage object stored in im
. Check the documentation for the available methods or check stackexchange for specific questions you might have.
You can also access the figure itself through ax.figure
. We will then save this to a file
ax.figure.savefig("colorful.tiff")
We can apply functions to an image to change it into another image. Any function that takes in an a 2D image and outputs another 2D image could be considered a filter
. You could write your own filters (currently there aren't so many implemented), and then you can use it with the apply_filter
command. Inside, you provide the filter function as first argument and then all the arguments and keyword arguments that the filter might need. Output is another image object.
Note that the function will always warn you that scales and units of the image may become incorrect, for example if the image is stretched or rebinned. By default the same pixelsize is kept, it is up to the user to evaluate whether the scale of the image should be updated.
gaussfilt = av.apply_filter(imf.gauss_filter, ks=10, sig=5, cval = av.data.mean())
medfilt = av.apply_filter(imf.med_filter)
butterfilt = av.apply_filter(imf.bw_filter)
WARNING:temmeta.data_io:The pixel size and unit of the original image was used. The scale may no longer be correct. Please verify and use set_scale. WARNING:temmeta.data_io:The pixel size and unit of the original image was used. The scale may no longer be correct. Please verify and use set_scale. WARNING:temmeta.data_io:The pixel size and unit of the original image was used. The scale may no longer be correct. Please verify and use set_scale.
We can use the same settings as before to plot some these images, for example the butterworth filtered image
butterfilt.plot(width = 10, sb_settings = scaleset, imshow_kwargs = imgset)
(<matplotlib.axes._axes.Axes at 0x149c19bb0>, <matplotlib.image.AxesImage at 0x149cab1c0>)
A few "filters" have been implemented as shorthand methods in the GeneralImage class already. These are currently:
Calling these methods, especially rebin, will automatically update the scale. You can add a cell and check with plotting that these functions worked
av_crop = av.crop(78, 72, 121, 171)
av_rebin = av.rebin(factor = 2)
av_linscale = av.linscale(min=2e4, max=5e4)
How do you know which values to take for minimum and maximum in linear scaling of the intensity? One way is using the plot_interactive
sliders to determine good thresholds. Another is to check the histogram of the image. This is easily done with plot_histogram
. With this you can check where the bulk of the intensity is located and adjust the scale accordingly.
av.plot_histogram()
(<AxesSubplot:xlabel='Pixel value (a.u.)', ylabel='Pixel count'>, <BarContainer object of 256 artists>)
You can also add, subtract and multiply images (element-wise) that are of the same size and scale. This is simply done with the +, - and * operators. However, it is important to be aware of over/underflow issues: if the image has the data type uint16
(most images in an emd file), it can store values between 0 and 65535. If two images are added and a pixel exceeds 65535, the value will loop back around to 0. To avoid issues like this, you can convert the data type of the image. Usually float is safe against overflow/underflow errors, but of course these take up much more RAM.
Changing the dtype is an in-place operation! If you change the dtype to something "lower" like uint8 which causes overflow errors, changing back to uint16 will not retrieve the original values!
av_linscale.change_dtype(float)
We demonstrate here a few random operations, which in truth don't make that much sense.
av_added = av_linscale + av_linscale
av_subtract = av_added - av_linscale
av_multiply = av_subtract*av_added
av_multiply.plot()
(<matplotlib.axes._axes.Axes at 0x14a77f160>, <matplotlib.image.AxesImage at 0x14abdb5e0>)
The metadata now becomes completely unreadable. This is because TEMMETA will try to keep all the metadata about all the images in nested structures. When you add images together, the resulting image is the child of two images, so TEMMETA will store all the info about both images in the child.
av_multiply.metadata
{ "data_axes": { "combos": [ { "navigation": null, "signal": [ "x", "y" ] } ], "x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" }, "y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" } }, "experiment_type": "modified", "parent_meta": [ { "data_axes": { "combos": [ { "navigation": null, "signal": [ "x", "y" ] } ], "x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" }, "y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" } }, "experiment_type": "modified", "parent_meta": [ { "data_axes": { "combos": [ { "navigation": null, "signal": [ "x", "y" ] } ], "x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" }, "y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" } }, "experiment_type": "modified", "parent_meta": [ { "data_axes": { "combos": [ { "navigation": null, "signal": [ "x", "y" ] } ], "x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" }, "y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" } }, "experiment_type": "modified", "parent_meta": { "data_axes": { "combos": [ { "navigation": null, "signal": [ "x", "y" ] } ], "x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" }, "y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" } }, "experiment_type": "modified", "parent_meta": { "acquisition": { "detector": { "collection_angles": { "begin": [ 25.41462224250735, "mrad" ], "end": [ 126.42261929381773, "mrad" ] }, "gain": [ 35.4, null ], "name": "HAADF", "offset": [ -2.6999999999999993, null ], "pixel_dwell_time": [ 1e-05, "s" ], "type": "ScanningDetector" }, "start_date_time": "2017-09-12T19:13:34" }, "data_axes": { "combos": [ { "navigation": [ "frame" ], "signal": [ "scan_x", "scan_y" ] }, { "navigation": [ "scan_x", "scan_y", "frame" ], "signal": null }, { "navigation": [ "time" ], "signal": null } ], "frame": { "axistype": "linear", "bins": 240, "offset": 0.0, "scale": 1.0, "unit": null }, "scan_x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 1.7610479458127432e-11, "unit": "m" }, "scan_y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 1.7610479458127432e-11, "unit": "m" }, "time": { "axistype": "function", "bins": 15728640, "function": "bin*1e-05+bin//256*0.0009000000000000002+bin//65536*9.99999999995449e-06", "unit": "s" } }, "experiment_type": "scan_image_series", "illumination": { "beam": { "convergence_angle": [ 17.027808627007598, "mrad" ], "mode": "convergent", "screen_current": [ 0.0, "A" ], "spot_index": 6 }, "source": { "acceleration_voltage": [ 300.0, "kV" ], "extraction_voltage": [ 3700.01220703125, "V" ] } }, "imaging": { "camera_length": [ 0.301, "m" ], "defocus": [ -2.4094582843801878e-08, "m" ], "mode": "diffraction", "stem_magnification": 20500000 }, "instrument": { "id": "3371" }, "meta": { "filename": "stem-edx.emd" }, "operator": {}, "optics": { "apertures": { "c2": { "order": "2", "position": { "x": [ 0.00442375, "m" ], "y": [ -0.000208125, "m" ] }, "shape": { "diameter": [ 4.9999999999999996e-05, "m" ] }, "type": "Circular" } }, "lenses": { "c1": { "current": [ -0.36628735065460205, null ] }, "c2": { "current": [ 0.213932603597641, null ] }, "c3": { "current": [ 0.35899609327316284, null ] }, "diffraction": { "current": [ 0.3123510479927063, null ] }, "intermediate": { "current": [ -0.5520547032356262, null ] }, "minicondenser": { "current": [ 0.18086114525794983, null ] }, "objective": { "current": [ 0.8855918049812317, null ] }, "projector1": { "current": [ -0.5004832148551941, null ] }, "projector2": { "current": [ -0.8716028928756714, null ] } } }, "process": "Extracted dataset Image/6fdbde41eecc4375b45cd86bd2be17c0 from EMD file stem-edx.emd", "sample": "", "software": { "name": "Velox", "version": "2.5.1" }, "stage": { "holder": "", "position": { "x": [ 0.00038286924, "m" ], "y": [ 0.000167710392, "m" ], "z": [ -0.00021097629000000001, "m" ] }, "tilt": { "alpha": [ -0.13979291438335673, "\u00b0" ], "beta": [ 0.2755002224273049, "\u00b0" ] } } }, "process": "Averaged all frames in stack" }, "process": "Scaled intensity between 20000.0 and 50000.0" }, { "data_axes": { "combos": [ { "navigation": null, "signal": [ "x", "y" ] } ], "x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" }, "y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" } }, "experiment_type": "modified", "parent_meta": { "data_axes": { "combos": [ { "navigation": null, "signal": [ "x", "y" ] } ], "x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" }, "y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" } }, "experiment_type": "modified", "parent_meta": { "acquisition": { "detector": { "collection_angles": { "begin": [ 25.41462224250735, "mrad" ], "end": [ 126.42261929381773, "mrad" ] }, "gain": [ 35.4, null ], "name": "HAADF", "offset": [ -2.6999999999999993, null ], "pixel_dwell_time": [ 1e-05, "s" ], "type": "ScanningDetector" }, "start_date_time": "2017-09-12T19:13:34" }, "data_axes": { "combos": [ { "navigation": [ "frame" ], "signal": [ "scan_x", "scan_y" ] }, { "navigation": [ "scan_x", "scan_y", "frame" ], "signal": null }, { "navigation": [ "time" ], "signal": null } ], "frame": { "axistype": "linear", "bins": 240, "offset": 0.0, "scale": 1.0, "unit": null }, "scan_x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 1.7610479458127432e-11, "unit": "m" }, "scan_y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 1.7610479458127432e-11, "unit": "m" }, "time": { "axistype": "function", "bins": 15728640, "function": "bin*1e-05+bin//256*0.0009000000000000002+bin//65536*9.99999999995449e-06", "unit": "s" } }, "experiment_type": "scan_image_series", "illumination": { "beam": { "convergence_angle": [ 17.027808627007598, "mrad" ], "mode": "convergent", "screen_current": [ 0.0, "A" ], "spot_index": 6 }, "source": { "acceleration_voltage": [ 300.0, "kV" ], "extraction_voltage": [ 3700.01220703125, "V" ] } }, "imaging": { "camera_length": [ 0.301, "m" ], "defocus": [ -2.4094582843801878e-08, "m" ], "mode": "diffraction", "stem_magnification": 20500000 }, "instrument": { "id": "3371" }, "meta": { "filename": "stem-edx.emd" }, "operator": {}, "optics": { "apertures": { "c2": { "order": "2", "position": { "x": [ 0.00442375, "m" ], "y": [ -0.000208125, "m" ] }, "shape": { "diameter": [ 4.9999999999999996e-05, "m" ] }, "type": "Circular" } }, "lenses": { "c1": { "current": [ -0.36628735065460205, null ] }, "c2": { "current": [ 0.213932603597641, null ] }, "c3": { "current": [ 0.35899609327316284, null ] }, "diffraction": { "current": [ 0.3123510479927063, null ] }, "intermediate": { "current": [ -0.5520547032356262, null ] }, "minicondenser": { "current": [ 0.18086114525794983, null ] }, "objective": { "current": [ 0.8855918049812317, null ] }, "projector1": { "current": [ -0.5004832148551941, null ] }, "projector2": { "current": [ -0.8716028928756714, null ] } } }, "process": "Extracted dataset Image/6fdbde41eecc4375b45cd86bd2be17c0 from EMD file stem-edx.emd", "sample": "", "software": { "name": "Velox", "version": "2.5.1" }, "stage": { "holder": "", "position": { "x": [ 0.00038286924, "m" ], "y": [ 0.000167710392, "m" ], "z": [ -0.00021097629000000001, "m" ] }, "tilt": { "alpha": [ -0.13979291438335673, "\u00b0" ], "beta": [ 0.2755002224273049, "\u00b0" ] } } }, "process": "Averaged all frames in stack" }, "process": "Scaled intensity between 20000.0 and 50000.0" } ], "process": "Added images element-wise" }, { "data_axes": { "combos": [ { "navigation": null, "signal": [ "x", "y" ] } ], "x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" }, "y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" } }, "experiment_type": "modified", "parent_meta": { "data_axes": { "combos": [ { "navigation": null, "signal": [ "x", "y" ] } ], "x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" }, "y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" } }, "experiment_type": "modified", "parent_meta": { "acquisition": { "detector": { "collection_angles": { "begin": [ 25.41462224250735, "mrad" ], "end": [ 126.42261929381773, "mrad" ] }, "gain": [ 35.4, null ], "name": "HAADF", "offset": [ -2.6999999999999993, null ], "pixel_dwell_time": [ 1e-05, "s" ], "type": "ScanningDetector" }, "start_date_time": "2017-09-12T19:13:34" }, "data_axes": { "combos": [ { "navigation": [ "frame" ], "signal": [ "scan_x", "scan_y" ] }, { "navigation": [ "scan_x", "scan_y", "frame" ], "signal": null }, { "navigation": [ "time" ], "signal": null } ], "frame": { "axistype": "linear", "bins": 240, "offset": 0.0, "scale": 1.0, "unit": null }, "scan_x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 1.7610479458127432e-11, "unit": "m" }, "scan_y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 1.7610479458127432e-11, "unit": "m" }, "time": { "axistype": "function", "bins": 15728640, "function": "bin*1e-05+bin//256*0.0009000000000000002+bin//65536*9.99999999995449e-06", "unit": "s" } }, "experiment_type": "scan_image_series", "illumination": { "beam": { "convergence_angle": [ 17.027808627007598, "mrad" ], "mode": "convergent", "screen_current": [ 0.0, "A" ], "spot_index": 6 }, "source": { "acceleration_voltage": [ 300.0, "kV" ], "extraction_voltage": [ 3700.01220703125, "V" ] } }, "imaging": { "camera_length": [ 0.301, "m" ], "defocus": [ -2.4094582843801878e-08, "m" ], "mode": "diffraction", "stem_magnification": 20500000 }, "instrument": { "id": "3371" }, "meta": { "filename": "stem-edx.emd" }, "operator": {}, "optics": { "apertures": { "c2": { "order": "2", "position": { "x": [ 0.00442375, "m" ], "y": [ -0.000208125, "m" ] }, "shape": { "diameter": [ 4.9999999999999996e-05, "m" ] }, "type": "Circular" } }, "lenses": { "c1": { "current": [ -0.36628735065460205, null ] }, "c2": { "current": [ 0.213932603597641, null ] }, "c3": { "current": [ 0.35899609327316284, null ] }, "diffraction": { "current": [ 0.3123510479927063, null ] }, "intermediate": { "current": [ -0.5520547032356262, null ] }, "minicondenser": { "current": [ 0.18086114525794983, null ] }, "objective": { "current": [ 0.8855918049812317, null ] }, "projector1": { "current": [ -0.5004832148551941, null ] }, "projector2": { "current": [ -0.8716028928756714, null ] } } }, "process": "Extracted dataset Image/6fdbde41eecc4375b45cd86bd2be17c0 from EMD file stem-edx.emd", "sample": "", "software": { "name": "Velox", "version": "2.5.1" }, "stage": { "holder": "", "position": { "x": [ 0.00038286924, "m" ], "y": [ 0.000167710392, "m" ], "z": [ -0.00021097629000000001, "m" ] }, "tilt": { "alpha": [ -0.13979291438335673, "\u00b0" ], "beta": [ 0.2755002224273049, "\u00b0" ] } } }, "process": "Averaged all frames in stack" }, "process": "Scaled intensity between 20000.0 and 50000.0" } ], "process": "Subtracted image 1 (right) from image 0 (left) element-wise" }, { "data_axes": { "combos": [ { "navigation": null, "signal": [ "x", "y" ] } ], "x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" }, "y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" } }, "experiment_type": "modified", "parent_meta": [ { "data_axes": { "combos": [ { "navigation": null, "signal": [ "x", "y" ] } ], "x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" }, "y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" } }, "experiment_type": "modified", "parent_meta": { "data_axes": { "combos": [ { "navigation": null, "signal": [ "x", "y" ] } ], "x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" }, "y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" } }, "experiment_type": "modified", "parent_meta": { "acquisition": { "detector": { "collection_angles": { "begin": [ 25.41462224250735, "mrad" ], "end": [ 126.42261929381773, "mrad" ] }, "gain": [ 35.4, null ], "name": "HAADF", "offset": [ -2.6999999999999993, null ], "pixel_dwell_time": [ 1e-05, "s" ], "type": "ScanningDetector" }, "start_date_time": "2017-09-12T19:13:34" }, "data_axes": { "combos": [ { "navigation": [ "frame" ], "signal": [ "scan_x", "scan_y" ] }, { "navigation": [ "scan_x", "scan_y", "frame" ], "signal": null }, { "navigation": [ "time" ], "signal": null } ], "frame": { "axistype": "linear", "bins": 240, "offset": 0.0, "scale": 1.0, "unit": null }, "scan_x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 1.7610479458127432e-11, "unit": "m" }, "scan_y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 1.7610479458127432e-11, "unit": "m" }, "time": { "axistype": "function", "bins": 15728640, "function": "bin*1e-05+bin//256*0.0009000000000000002+bin//65536*9.99999999995449e-06", "unit": "s" } }, "experiment_type": "scan_image_series", "illumination": { "beam": { "convergence_angle": [ 17.027808627007598, "mrad" ], "mode": "convergent", "screen_current": [ 0.0, "A" ], "spot_index": 6 }, "source": { "acceleration_voltage": [ 300.0, "kV" ], "extraction_voltage": [ 3700.01220703125, "V" ] } }, "imaging": { "camera_length": [ 0.301, "m" ], "defocus": [ -2.4094582843801878e-08, "m" ], "mode": "diffraction", "stem_magnification": 20500000 }, "instrument": { "id": "3371" }, "meta": { "filename": "stem-edx.emd" }, "operator": {}, "optics": { "apertures": { "c2": { "order": "2", "position": { "x": [ 0.00442375, "m" ], "y": [ -0.000208125, "m" ] }, "shape": { "diameter": [ 4.9999999999999996e-05, "m" ] }, "type": "Circular" } }, "lenses": { "c1": { "current": [ -0.36628735065460205, null ] }, "c2": { "current": [ 0.213932603597641, null ] }, "c3": { "current": [ 0.35899609327316284, null ] }, "diffraction": { "current": [ 0.3123510479927063, null ] }, "intermediate": { "current": [ -0.5520547032356262, null ] }, "minicondenser": { "current": [ 0.18086114525794983, null ] }, "objective": { "current": [ 0.8855918049812317, null ] }, "projector1": { "current": [ -0.5004832148551941, null ] }, "projector2": { "current": [ -0.8716028928756714, null ] } } }, "process": "Extracted dataset Image/6fdbde41eecc4375b45cd86bd2be17c0 from EMD file stem-edx.emd", "sample": "", "software": { "name": "Velox", "version": "2.5.1" }, "stage": { "holder": "", "position": { "x": [ 0.00038286924, "m" ], "y": [ 0.000167710392, "m" ], "z": [ -0.00021097629000000001, "m" ] }, "tilt": { "alpha": [ -0.13979291438335673, "\u00b0" ], "beta": [ 0.2755002224273049, "\u00b0" ] } } }, "process": "Averaged all frames in stack" }, "process": "Scaled intensity between 20000.0 and 50000.0" }, { "data_axes": { "combos": [ { "navigation": null, "signal": [ "x", "y" ] } ], "x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" }, "y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" } }, "experiment_type": "modified", "parent_meta": { "data_axes": { "combos": [ { "navigation": null, "signal": [ "x", "y" ] } ], "x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" }, "y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 0.017610479458127434, "unit": "nm" } }, "experiment_type": "modified", "parent_meta": { "acquisition": { "detector": { "collection_angles": { "begin": [ 25.41462224250735, "mrad" ], "end": [ 126.42261929381773, "mrad" ] }, "gain": [ 35.4, null ], "name": "HAADF", "offset": [ -2.6999999999999993, null ], "pixel_dwell_time": [ 1e-05, "s" ], "type": "ScanningDetector" }, "start_date_time": "2017-09-12T19:13:34" }, "data_axes": { "combos": [ { "navigation": [ "frame" ], "signal": [ "scan_x", "scan_y" ] }, { "navigation": [ "scan_x", "scan_y", "frame" ], "signal": null }, { "navigation": [ "time" ], "signal": null } ], "frame": { "axistype": "linear", "bins": 240, "offset": 0.0, "scale": 1.0, "unit": null }, "scan_x": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 1.7610479458127432e-11, "unit": "m" }, "scan_y": { "axistype": "linear", "bins": 256, "offset": 0.0, "scale": 1.7610479458127432e-11, "unit": "m" }, "time": { "axistype": "function", "bins": 15728640, "function": "bin*1e-05+bin//256*0.0009000000000000002+bin//65536*9.99999999995449e-06", "unit": "s" } }, "experiment_type": "scan_image_series", "illumination": { "beam": { "convergence_angle": [ 17.027808627007598, "mrad" ], "mode": "convergent", "screen_current": [ 0.0, "A" ], "spot_index": 6 }, "source": { "acceleration_voltage": [ 300.0, "kV" ], "extraction_voltage": [ 3700.01220703125, "V" ] } }, "imaging": { "camera_length": [ 0.301, "m" ], "defocus": [ -2.4094582843801878e-08, "m" ], "mode": "diffraction", "stem_magnification": 20500000 }, "instrument": { "id": "3371" }, "meta": { "filename": "stem-edx.emd" }, "operator": {}, "optics": { "apertures": { "c2": { "order": "2", "position": { "x": [ 0.00442375, "m" ], "y": [ -0.000208125, "m" ] }, "shape": { "diameter": [ 4.9999999999999996e-05, "m" ] }, "type": "Circular" } }, "lenses": { "c1": { "current": [ -0.36628735065460205, null ] }, "c2": { "current": [ 0.213932603597641, null ] }, "c3": { "current": [ 0.35899609327316284, null ] }, "diffraction": { "current": [ 0.3123510479927063, null ] }, "intermediate": { "current": [ -0.5520547032356262, null ] }, "minicondenser": { "current": [ 0.18086114525794983, null ] }, "objective": { "current": [ 0.8855918049812317, null ] }, "projector1": { "current": [ -0.5004832148551941, null ] }, "projector2": { "current": [ -0.8716028928756714, null ] } } }, "process": "Extracted dataset Image/6fdbde41eecc4375b45cd86bd2be17c0 from EMD file stem-edx.emd", "sample": "", "software": { "name": "Velox", "version": "2.5.1" }, "stage": { "holder": "", "position": { "x": [ 0.00038286924, "m" ], "y": [ 0.000167710392, "m" ], "z": [ -0.00021097629000000001, "m" ] }, "tilt": { "alpha": [ -0.13979291438335673, "\u00b0" ], "beta": [ 0.2755002224273049, "\u00b0" ] } } }, "process": "Averaged all frames in stack" }, "process": "Scaled intensity between 20000.0 and 50000.0" } ], "process": "Added images element-wise" } ], "process": "Multiplied images element-wise" }
The best way to visualize what has happened to the image is using the history
attribute. When the dataset derives from multiple images you will see a tree structure, with the most recent operations at the top, including all operations on all parents in different branches. You can verify that this tree corresponds to our operations higher up
av_multiply.history
* Multiplied images element-wise |\ |* Added images element-wise ||\ ||* Scaled intensity between 20000.0 and 50000.0 ||| ||* Averaged all frames in stack ||| ||* Extracted dataset Image/6fdbde41eecc4375b45cd86bd2be17c0 from EMD file stem-edx.emd || |* Scaled intensity between 20000.0 and 50000.0 || |* Averaged all frames in stack || |* Extracted dataset Image/6fdbde41eecc4375b45cd86bd2be17c0 from EMD file stem-edx.emd | * Subtracted image 1 (right) from image 0 (left) element-wise |\ |* Scaled intensity between 20000.0 and 50000.0 || |* Averaged all frames in stack || |* Extracted dataset Image/6fdbde41eecc4375b45cd86bd2be17c0 from EMD file stem-edx.emd | * Added images element-wise |\ |* Scaled intensity between 20000.0 and 50000.0 || |* Averaged all frames in stack || |* Extracted dataset Image/6fdbde41eecc4375b45cd86bd2be17c0 from EMD file stem-edx.emd | * Scaled intensity between 20000.0 and 50000.0 | * Averaged all frames in stack | * Extracted dataset Image/6fdbde41eecc4375b45cd86bd2be17c0 from EMD file stem-edx.emd
Finally, we can also export our images to png, tiff, and hyperspy. With the save
command, the image is by default exported as an 8-bit png. However, you can add options min
, max
and dtype
to change this behavior, checkout help(av.save)
for this. There is no scalebar applied.
For exporting as tiff (by default with the same datatype as stored in the original dataset), we use the plot command as seen before, and export the plot by providing the filename
argument.
For creating a hyperspy dataset, we follow the same example as before.
av.save("average.png")
av.to_hspy("average")
<Signal2D, title: , dimensions: (|256, 256)>
We can import files again with import_file_to_image
. The function will also look for a file ending in x_meta.json
in the same folder, with x
the name of the image file. Normally this is automatically written out when you plot or save an image.
av_import = dio.import_file_to_image("average.png")
We can calculate the fast fourrier transform of images as long as they are square. You can actually calculate an FFT of an image of any size, but a non-square one is not so intuitive to interpret, so we don't allow it in TEMMETA. We can do a number of things on the FFT, mainly extracting fourrier components by masking and calculating the inverse fourrier transform. We can also calculate strain with GPA. These operations are demonstrated here on an HRTEM image dataset. We first import the emd file and extract the dataset
emd3 = dio.EMDFile(emd3_path)
emd3.print_simple_structure()
------------------------------ Image datasets ------------------------------ Dataset number: 0, UUID: 7c3929705ae04354b656fb2f929a1cfa, Shape: 4096x4096, Frames: 1, Data type: int16, Min:179, Max:1966
hrtemimg = emd3.get_dataset("Image", "7c3929705ae04354b656fb2f929a1cfa")
hrtemimg.plot()
(<matplotlib.axes._axes.Axes at 0x14abed5b0>, <matplotlib.image.AxesImage at 0x14b715670>)
Scales are again in meters so we set it to nm
hrtemimg.set_scale(hrtemimg.pixelsize*1e9, "nm")
We access the fft with
fft = hrtemimg.fft
For quick visualisation of the FFT and for things like coordinate selection for creating masks, there is the plot function. Note however that an FFT is complex valued so what is actually plot is a power spectrum on a log scale. The function will warn you of this. For optimizing the image of the FFT in case you want to write this out, make an image from the real or imaginary component, or from the argument or phase.
ax, im = fft.plot()
im.set_clim(0.7, 1)
WARNING:temmeta.data_io:An FFT is complex-valued and can not be visualized as an image. Actually plot here is a smoothed power spectrum, logarithmically scaled. If the FFT in point (x, y) is a+bi, then the image intensity there is log10(a^2+b^2)
For example, I might want this same image but without a scale bar, with some modified intesity scaling, and with a different colormap. Below is an example. In the first cell we extract the power spectrum (log scaled), then we apply a median filter, then we plot with a few additional parameters
ps = fft.power_spectrum
ps.apply_filter(imf.med_filter, inplace=True)
WARNING:temmeta.data_io:The pixel size and unit of the original image was used. The scale may no longer be correct. Please verify and use set_scale. WARNING:temmeta.data_io:Data in object will be overwritten
msize = 800 # some variable for the window size
ax, im = ps.plot(scale_bar=False, width=20, imshow_kwargs={"cmap": "hot", "vmin": 100, "vmax": 150})
ax.set_xlim(ps.width//2-msize//2, ps.width//2+msize//2)
# note that to keep the same orientation as the image, we have to go from positive to negative on the y-axis
# This is because when plotting images, the y-axis points from the top left corner downwards.
ax.set_ylim(ps.height//2+msize//2, ps.height//2-msize//2)
(2448.0, 1648.0)
Saving this figure
ax.figure.savefig("fft.png")
So much for visualizing and plotting the fft.
Now we will try to do some more operations on the FFT stored in fft
. You can measure distances and angles on the FFT like in images beforehand. There is a small difference, in that if you provide only one point, it will automatically take the central spot as the second point. You can see this in the following cell
print(fft.measure(2151, 2140))
print(fft.measure(fft.width/2, fft.height/2, 2151, 2140))
((4.909470530547582, '1/nm'), (41.78779023989079, '°')) ((4.884501060243416, '1/nm'), (41.77135213036438, '°'))
We can get an ifft from the FFT via ifft
. If we get the ifft from our fft without processing we retrieve our original image.
Note: The IFFT is also a complex valued array in general. If the FFT is symmetric over 0 then the IFFT will be real-valued except for small rounding errors. Still, what we actually return as an FFT is the modulus. So if you invert a non-centrosymmetric FFT, the ifft is complex valued. But in this implementation, we always take the modulus to work with and plot a real image.
ifft = fft.ifft
ifft.plot()
(<matplotlib.axes._axes.Axes at 0x14ab08790>, <matplotlib.image.AxesImage at 0x14ace85e0>)
You should know however that it is not completely the same. Due to the calculations, the values are converted to floating point numbers with decimals. You will see this when you check the dtype
of the data. Due to various rounding errors there will be a slight but negligible difference between the values in each pixel.
ifft.data.dtype
dtype('float64')
The metadata will keep track of these operations
ifft.history
* Inverse Fourrier Transform | * Fast-fourrier-transform | * Extracted dataset Image/7c3929705ae04354b656fb2f929a1cfa from EMD file hrtem.emd
We will now create some masks for the FFT. This allows us to select/exclude certain spatial frequencies. We can then see the effect on the image with an inverse fourrier transform. A mask is nothing but an array of the same size as the FFT with values between 0 and 1. Multiplying the mask element-wise with the FFT will act as a frequency filter. Any grayscale image of the same size as the FFT can act as a mask. With TEMMETA we provide a wrapper class Mask, which includes a few handy functions for adding elements to it.
In TEMMETA, a Mask is initialized as either a positive or a negative mask. A positive mask is where all the elements are intialized as 1. If no features are added, everything passes the filter. Features added are "black", so you add features to exclude frequencies from the FFT. A negative mask is initilized at 0 and nothing is passed. Features added are "white", so these are made to allow some frequencies through. For high resolution microscopy, negative masks with a few "holes" are the most common, so this is the default behavior. If you want to change the default behavior use the argument negative=False
. You can also always invert a mask with the method invert()
.
We can instantiate a Mask
object directly from the FFT.
m1 = fft.create_mask()
We then use the methods of Mask
to add features to it. Currently implemented are (check their help for details):
mirrored=False
.mirrored
can be used.mirrored
can be used. The center is not blanked.In the case of a negative mask (all is 0), you "poke holes" (where value is 1) in it of various shapes. In the case of a positive mask (all is 1), you "paint features" (where value is 0) over it. You can also paint features in negative masks or poke holes in positive masks, for example to poke a hole in a painted feature. All these things are illustrated below.
# some partially overlapping circles which are mirrored
m1.add_circle(2107, 1923, 25, edge_blur=5)
m1.add_circle(2150, 1923, 25, edge_blur=5)
m1.add_circle(2150, 1970, 50, edge_blur=15)
# paint a small circle on a previous hole
m1.add_circle(2150, 1923, 5, edge_blur=0, negative=True)
# a square which is not mirrored
m1.add_square(2150, 2142, 20, edge_blur=2, mirrored=False)
# a pass-band
m1.add_band(20, 50, edge_blur=2)
# a 2D array
m1.add_array2D(2335, 2124, 1981, 1614, 10, edge_blur=0)
Also the history of masks is kept
m1.history
* Added positive 2D array. Radius: 10. Positions: (2335,2124), (1981,1614). Edge_blur: 0. | * Added positive pass-band. Radii: 20, 50. Edge_blur: 2. | * Added positive square. Side-length: 40. Position: (2150, 2142). Mirrored: False. Edge_blur: 2. | * Added negative circle. Radius: 5. Position: (2150, 1923). Mirrored: True. Edge_blur: 0. | * Added positive circle. Radius: 50. Position: (2150, 1970). Mirrored: True. Edge_blur: 15. | * Added positive circle. Radius: 25. Position: (2150, 1923). Mirrored: True. Edge_blur: 5. | * Added positive circle. Radius: 25. Position: (2107, 1923). Mirrored: True. Edge_blur: 5. | * Initialized negative mask of size 4096x4096
You can visualize the mask directly with plot
. Use the zoom button to inspect the mask in greater detail.
m1.plot()
(<AxesSubplot:>, <matplotlib.image.AxesImage at 0x14ac8d820>)
You can invert the mask to turn it from positive to negative and vice versa
m1.invert()
m1.plot()
(<AxesSubplot:>, <matplotlib.image.AxesImage at 0x14b631d60>)
I hope with this example you see you have a lot of options to designing an appropriate FFT mask. Of course this example is a bit nonsense, so we will create one relevant to our fft with a 3 selected reflections. I measured the (x,y) coordinates on a plot of the fft, then enter them into the function below. I also add the argument refine=True
, which will optimize the positioning of the feature so it is centered on the maximum intensity pixel.
m2 = fft.create_mask()
m2.add_circle(2204, 2012, 15, refine=True)
m2.add_circle(2150, 2140, 15, refine=True)
m2.add_circle(2108, 1916, 15, refine=True)
#m2.add_array2D(2204, 2012, 2147, 2140, 10)
WARNING:temmeta.data_io:Refined coordinates: (2204, 2013) WARNING:temmeta.data_io:Refined coordinates: (1892, 2083) WARNING:temmeta.data_io:Refined coordinates: (2148, 2140) WARNING:temmeta.data_io:Refined coordinates: (1948, 1956) WARNING:temmeta.data_io:Refined coordinates: (2105, 1922) WARNING:temmeta.data_io:Refined coordinates: (1991, 2174)
I can check the mask on top of the fft with check_mask
or check_masked
. The first plots the mask as a partially transparent image over the power spectrum. The second plots the power spectrum if the fft is masked. With the values that are returned from the plot function I can also optimize the plot a bit.
ax, imfft, immask = fft.check_mask(m2)
imfft.set_clim(100, 150)
When we are happy with our mask we can apply it to the FFT and then take the inverse Fourrier transform
fft_masked = fft.apply_mask(m2)
ifft_masked = fft_masked.ifft
ifft_masked
is again an image which we can use for the same kinds of operations as before
ifft_masked.plot()
(<matplotlib.axes._axes.Axes at 0x14b659af0>, <matplotlib.image.AxesImage at 0x14bd318e0>)
Once again, the metadata has kept a record of the operations. You will see a similar branched structure as when adding/subtracting multiple images.
ifft_masked.history
* Inverse Fourrier Transform | * Applied mask (multiplied) |\ |* Added positive circle. Radius: 15. Position: (2105, 1922). Mirrored: True. Edge_blur: 3. || |* Added positive circle. Radius: 15. Position: (2148, 2140). Mirrored: True. Edge_blur: 3. || |* Added positive circle. Radius: 15. Position: (2204, 2013). Mirrored: True. Edge_blur: 3. || |* Initialized negative mask of size 4096x4096 | * Fast-fourrier-transform | * Extracted dataset Image/7c3929705ae04354b656fb2f929a1cfa from EMD file hrtem.emd
mimg = m2.to_image()
mimg.save("mask.png", dtype=np.uint8)
I also implemented a very basic version of geometric phase analysis (GPA) (Hytch et al., 1998), which calculates strain maps in HRTEM images. The application is pretty basic: use the method calc_strains
on the image. You must provide as arguments the coordinates of two spots in the FFT, from which phase images are internally calculated and processed. For the other arguments, check the documentation with help
. For a large image as this one the calculation can take quite some time, half a minute to a minute or two.
The method returns 4 image objects representing different strain components: epsilon_xx, epsilon_yy, epsilon_xy, the components from the symmetric strain matrix, and omega_xy, the off-diagonal element of the antisymmetric term representing rigid rotation. These can be plot and worked with as any other image, see previously. Note also that since this method is based on derivatives of the phase, the output images are 1 pixel less wide and high than the original.
An additional argument that can be used is angle
, which describes the strain tensor in another coordinate system. If your desired x-axis lies 30 degrees with the horizontal for example, angle=30
and epsilon_xx will represent the tensile/compressive strain in this direction.
The method must still be verified against other implementations and better datasets, don't trust the result yet for your publications
Here a demonstration of how it's supposed to work. we will use a subset of hrtemimg from before. We use crop and to keep it square and a power of two, we add the argument nearest_power=True.
imgcropped = hrtemimg.crop(1809,2077,3300,3381, nearest_power=True)
imgcropped.plot()
(<matplotlib.axes._axes.Axes at 0x14bd462e0>, <matplotlib.image.AxesImage at 0x14b698310>)
We now check the FFT and find the coordinates of spots we want to select
fft_crop = imgcropped.fft
ax, im = fft_crop.plot()
im.set_clim(0.7, 1)
WARNING:temmeta.data_io:An FFT is complex-valued and can not be visualized as an image. Actually plot here is a smoothed power spectrum, logarithmically scaled. If the FFT in point (x, y) is a+bi, then the image intensity there is log10(a^2+b^2)
xx1, yy1 = 536, 536
xx2, yy2 = 550, 503
exx, eyy, exy, oxy = imgcropped.calc_strains(xx1, yy1, xx2, yy2, r=8, edge_blur=5)
exx.history
* Calculated a GPA map. Reflections: (536, 536), (550, 503). Radius: 8. Edge_blur: 5. Angle: 0 °. Epsilon_xx component. | * Cropped between x=1809-2833 and y=2077-3101 | * Extracted dataset Image/7c3929705ae04354b656fb2f929a1cfa from EMD file hrtem.emd
To make a somewhat more complex plot and put all images on the same figure, I access the data inside these GeneralImage objects directly and use matplotlib. In this way you can check how changing the parameters in the cell above changes the output. The result in general looks pretty crappy, but then again this isn't the best dataset.
fig, axes = plt.subplots(2, 2, figsize = (10,10))
(ax1, ax2), (ax3, ax4) = axes
kwargs = {"cmap": "hot", "vmin":-0.1, "vmax":0.1}
# exx
ax1.imshow(exx.data, **kwargs)
ax1.set_title(r"$\epsilon_{xx}$")
ax1.set_axis_off()
# eyy
ax2.imshow(eyy.data, **kwargs)
ax2.set_title(r"$\epsilon_{yy}$")
ax2.set_axis_off()
# exy
ax3.imshow(exy.data, **kwargs)
ax3.set_title(r"$\epsilon_{xy}$")
ax3.set_axis_off()
# oxy
ax4.imshow(oxy.data, **kwargs)
ax4.set_title(r"$\omega_{xy}$")
ax4.set_axis_off()
You can calculate the positions of atomic columns in HR STEM images using the find_peaks
function in the fitPeaks
module. To use it you must import it. I want to thank C. Liebscher for contributing this piece.
from temmeta import fitPeaks as fp
Finding can be a bit tricky so there are quite a few parameters and processing steps which happen in the background. Check the help of the function for more details. Currently these functions are not integrated in TEMMETA objects, so you have to pass the raw data of the image to the function. This will probably change at some point if there is demand for it.
pks = fp.find_peaks(av.data, gauss_args=(3, 0.5), min_dist=3, min_int=0.05)
pks
x | y | intensity | |
---|---|---|---|
0 | 17.0 | 26.0 | 0.098586 |
1 | 69.0 | 11.0 | 0.137608 |
2 | 33.0 | 21.0 | 0.150329 |
3 | 154.0 | 170.0 | 0.153108 |
4 | 159.0 | 185.0 | 0.154511 |
... | ... | ... | ... |
468 | 130.0 | 49.0 | 0.783765 |
469 | 57.0 | 163.0 | 0.795588 |
470 | 34.0 | 78.0 | 0.804434 |
471 | 73.0 | 185.0 | 0.808195 |
472 | 33.0 | 246.0 | 0.889151 |
473 rows × 3 columns
As you can see, the peaks are stored in an x, y intensity table. You can plot them quickly over the image as seen below. At some point there may come a wrapper function.
fig, ax = plt.subplots()
ax.imshow(av.data)
ax.scatter(pks["x"], pks["y"], color = "red")
<matplotlib.collections.PathCollection at 0x15911a6a0>
Once you have peaks, it is also possible to fit gaussians to the atomic columns and refine the peak position to sub-pixel accuracy with gaussian_peak_fit
. This process can be quite slow so beware. You must at least pass in the data and the starting peaks. A table will be returned with the fitted gaussians.
rpks = fp.gaussian_peak_fit(av.data, pks)
rpks
x | y | intensity | sigma | x refined | y refined | background | integrated intensity | |
---|---|---|---|---|---|---|---|---|
0 | 17.0 | 26.0 | 0.098586 | 3.200000 | 17.298447 | 26.186821 | 1114.853693 | 1485084.0 |
1 | 69.0 | 11.0 | 0.137608 | 3.200000 | 69.399562 | 11.359356 | 1403.339826 | 1575719.0 |
2 | 33.0 | 21.0 | 0.150329 | 3.200000 | 33.083573 | 21.354752 | -882.422611 | 1413470.0 |
3 | 154.0 | 170.0 | 0.153108 | 3.200000 | 154.013457 | 170.208414 | 18002.808806 | 2044844.0 |
4 | 159.0 | 185.0 | 0.154511 | 3.313547 | 158.415839 | 185.080843 | 13644.143556 | 2195072.0 |
... | ... | ... | ... | ... | ... | ... | ... | ... |
468 | 130.0 | 49.0 | 0.783765 | 3.200000 | 129.590476 | 49.136299 | 6525.704960 | 1888741.0 |
469 | 57.0 | 163.0 | 0.795588 | 5.168145 | 56.424296 | 163.556632 | -7140.133746 | 1614178.0 |
470 | 34.0 | 78.0 | 0.804434 | 3.200000 | 33.495185 | 77.774333 | 1999.091973 | 1824444.0 |
471 | 73.0 | 185.0 | 0.808195 | 3.528082 | 72.952925 | 185.255291 | 15762.040065 | 1595106.0 |
472 | 33.0 | 246.0 | 0.889151 | 3.200000 | 33.342749 | 245.537548 | 15730.212626 | 2131117.0 |
473 rows × 8 columns
fig, ax = plt.subplots()
ax.imshow(av.data)
ax.scatter(rpks["x refined"], rpks["y refined"], color = "red")
<matplotlib.collections.PathCollection at 0x159032970>
At the moment, the peak finding algorithm does not work well for FFTs or diffraction patterns!
Image stacks arrise mainly due to multiple frames being captured over time. We already showed a little bit of what one can do with a stack. The bulk of the analysis in principle happens in individual images, so you would convert frames or average multiple frames into one image and then work with it. But some operations you might like to do on all the images in the stack. Currently, the options are a bit limited but may expand as more of people's needs come to light.
Let's revisit the stack from before stored in stack
. We can either access individual frames with get_frame
or we can create new stacks with select_frames
or exclude_frames
frame203 = stack.get_frame(203)
frames_select = [1, 3, 5]
stack_select = stack.select_frames(frames_select)
stack.exclude_frames(frames_select, inplace=True) # to save on memory we do this inplace
If you select only one frame the returned object is an image which you can process in the same ways as before. If you select frames, returned is a new stack which you can further process.
To select areas you can use the crop
command, with the same behavior as for images, only that here it is applied to all frames. You could again use the nearest_power
argument to make sure your selection is a square of size $2^n \times 2^n$
stack.crop(50, 50, 120, 75, inplace=True)
The basic filters linscale
and rebin
also work on an image stack. There is also apply_filter
which just loops over all the frames and applies the filter. Some of these operations can be quite slow. There is an option to use multiprocessing to speed up the calculations a bit - don't expect magic though
stack = emd1.get_dataset("Image", "6fdbde41eecc4375b45cd86bd2be17c0")
stack.linscale(1e4, 4e4, inplace=True)
stack.rebin(1.5, multiprocessing=True, inplace=True)
stack = emd1.get_dataset("Image", "6fdbde41eecc4375b45cd86bd2be17c0")
# here we demonstrate the minor differences between using multiprocessing and not. (key word) args are just
# passed directly to the function.
from time import time
start = time()
stack_gausfilter = stack.apply_filter(imf.gauss_filter, multiprocessing=False, sig=5)
stop = time()
print(f"{(stop-start)} seconds without multiprocessing")
start = time()
stack_gausfilter = stack.apply_filter(imf.gauss_filter, multiprocessing=True, sig=5)
stop = time()
print(f"{(stop-start)} seconds with multiprocessing")
WARNING:temmeta.data_io:The pixel size and unit of the original images were used. The scale may no longer be correct. Please verify and use set_scale.
1.3539018630981445 seconds without multiprocessing
WARNING:temmeta.data_io:The pixel size and unit of the original images were used. The scale may no longer be correct. Please verify and use set_scale.
1.1902351379394531 seconds with multiprocessing
stack_gausfilter.plot_interactive()
Aligning features in different images is a process called image registration. If it is assumed the only difference between the images are shifts are rotations, aligning the images is called rigid registration. If there are deformations between the images, aligning features in the images requires destoring the images and is a very complex process called non-rigid-registration.
A very basic rigid registration algorithm (only correcting translations) is implemented in TEMMETA and can be applied on a stack using calc_misalignment
and align_frames
. It works by selecting a small square in the images and calculating the cross-correlation between this square and a search area in consecutive images or the first image. The maximum in the cross-correlation may be shifted from the center and this shift is applied on the frames. No deformations in individual frames can be corrected this way!
If you require non rigid registration (for example for strain measurements) you may want to look at match-series
implemented by Benjamin Berkels. Unfortunately since this code is written in C++ I can't easily integrate it into TEMMETA, so you will have to do a bit more work. I have written a few scripts which should make it a bit easier to extract the information inside an EMD file, run the non-rigid registration on it. Once the non-rigid registration has been applied you can import the images/spectra again and continue analyzing them with TEMMETA or another tool of your choice like Hyperspy
We will re-use the same image stack from before, although these are already pretty well aligned. Aligning periodic images, especially noisy ones like these, can be pretty tough because the cross-correlation can give multiple maxima. First we define a rectangle that we want to "follow" - this is given in red. Don't use the entire image: the cross correlations will take forever to calculate and it will almost certainly not find the optimum. Then we also define a "border thickness", a region in which we will look for the red region.
border = 8
rectangle = (75, 30, 105, 65)
rectangle2 = [75-border, 30-border, 105+border, 65+border]
fig, ax = plt.subplots()
ax.imshow(stack.data[0])
rx = [rectangle[0], rectangle[2], rectangle[2], rectangle[0], rectangle[0]]
ry = [rectangle[1], rectangle[1], rectangle[3], rectangle[3], rectangle[1]]
rx2 = [rectangle2[0], rectangle2[2], rectangle2[2], rectangle2[0], rectangle2[0]]
ry2 = [rectangle2[1], rectangle2[1], rectangle2[3], rectangle2[3], rectangle2[1]]
ax.plot(rx, ry, color="red")
ax.plot(rx2, ry2, color=(0, 1, 0))
[<matplotlib.lines.Line2D at 0x149c08130>]
We use calc_misalignment
to calculate the x and y shifts required to match the red region in all images. We pass the rectangle and border as arguments. In addition, there are 2 ways the method can run: all frames are compared to the first frame (cumulative=False (default)), or all frames are compared to the frame before it and the shifts are added consecutively (cumulative=True). In principle, comparing to the frame just before is better because these images will be most similar. However, errors will accumulate. Since we only have small shifts, we will set this to false. reg
is for regularisation: maxima further away from the center are penalized and a higher regularization term increases the importance of distance from the center. If reg
=0, there is no distance penalisation, which we use here. With debug=True
, we will also return the cross correlations as an image-stack which we can explore for debugging purposes.
For the entire stack, the calculation can take quite a while! The larger the area of interest and the larger the stack, the longer it takes.
(shifts, ccx) = stack.calc_misalignment(rectangle=rectangle, border=border, cumulative=False,
debug=True, reg=0)
You can visualize the cross correlations with plot_interactive
. In this case the brightest pixel is always more or less in the middle and shifts are almost 0. You can visualize it best by sliding the min slider higher.
ccx.plot_interactive()
We can apply the deformations with align_frames
, simply passing in the shifts calculated before. If no shifts are passed in, it will attempt to calculate the shifts in the same command, however I don't recommend this as you don't control the process as well. I have just cropped the image to the size of the search window for visualisation. You can see that the images still don't feel very aligned; this is mainly due to non-linear distortions which can not be corrected this way.
stack.align_frames(shifts = shifts, inplace=True)
stack.crop(*rectangle2, inplace=True)
stack.plot_interactive()
Image registration is an active field of research and I strongly welcome additions/improvements to the stack alignment methods.
There are currently three options for exporting images in a stack: export frames directly using export_frames
(fastest), plot frames with a scalebar and various options using plot_export_frames
(much slower), and converting and saving to a hyperspy dataset with to_hspy
First we will define some folder where we will export the images to. The export commands will create the folder if it doesn't exist. This does not happen with to_hspy
export_folder = "./Test_stack"
export_frames
is the easiest, just directly saving the raw image to a png with the same dtype as the stored image (unless otherwise specified). It will give you an estimate for the time it took to execute the command
stack.export_frames(export_folder+"/simple/", "Frame")
'export_frames' 0.69 s
For very large stacks, it's possible that multithreading can speed up the export a bit. By default it is turned off but you can turn it on with the option multithreading
. counter
indicates how many counter digits are appended to the frame. name
indicates the filename prefix (the part before the number)
stack.export_frames(export_folder, name = "boo", multithreading=True, counter = 4)
'export_frames' 0.35 s
If you don't want to export all frames, you can pass a list of frames to export
stack.export_frames(export_folder, name = "boo", counter = 4, frames = [2, 5, 6])
'export_frames' 0.05 s
If you want to export all the frames with a scalebar or a special color scheme, each one has to be plot and processed. The command works in a similar way but you can add additional arguments passed to plot_image
, such as settings for the image plot and settings for the scalebar. This command will take a lot longer than just exporting the raw images. For large stacks, multiprocessing
may give a slight time gain. For small arrays this will not yield time gains, it may even take longer.
imskw = {"cmap": "jet", "vmin": 2.5e4, "vmax": 5e4}
scaleset = {"location": 'lower right',
"color": "white",
"length_fraction": 0.25,
"frameon": True,
"box_color": "black",
"font_properties": {"size": 14, "weight": "bold"}}
# frames = [10, 20, 30, 60, 66, 88]
frames = None
stack.plot_export_frames(export_folder, name="plotted", multiprocessing=True, counter = 4,
frames = frames, width=10, imshow_kwargs=imskw, sb_settings = scaleset)
'plot_export_frames' 8.63 s
Finally, you can export your data to a hyperspy dataset which you can continue to process with hyperspy
stack_hspy = stack.to_hspy(filepath = export_folder+"/stack_hs")
You might already have a folder with images and wish to re-import them as a stack. Or you might want to stack a list of individual images. These are possible with the dio.import_files_to_stack
and dio.images_to_stack
. Note that this operation only makes sense if the images are related via the same scale.
With images_to_stack
, the scale of the first image is adopted as the scale and the metadata is copied.
If the images are read from a folder, with dio.import_files_to_stack
, the files must follow a specific naming convention in order to re-create the stack in the correct order: somearbitrarynamewithoutspaces_000x.png or somearbitrarynamewithoutspaces_x.png with x a number. All images should have the same number of counter digits. It is recommended to not have any other files in the folder except a metadata.json
file.
stack_import = dio.import_files_to_stack(export_folder+"/simple/")
stack_import.plot_interactive()
frames = [stack.get_frame(i) for i in [12, 35, 67, 89]]
restack = dio.images_to_stack(frames)
restack.plot_interactive()
TEMMETA can at present only read SpectrumStream data from EMD files, that is spectrum maps taken over multiple frames. By adding together the frames, one can obtain a spectrum map dataset. Such a dataset is what you will want to work with most for analysis.
Keeping the data in the stream format is convenient for one main reason: data from individual frames can be separated, adding a time component to the dataset. Hence, one can select/exclude frames, or as in the case of image stacks one can align each individual frame before adding together all the data. Rigid and non-rigid registration corrections, calculated from an image stack, can be applied to a spectrum stream, before adding all the frames together.
The downside is that Spectrum streams are extremely sparse: most scan positions do not yield any data and so are just 0; it is only by collecting multiple frames that the dataset fills up. Fitting an edx map acquired over an hour for hundreds of frames entirely into memory is very wastefull. Therefore SpectrumStream
stores its data in a sparse scipy.spmatrix
matrix form. This matrix has the shape (x*y*frames, channels), with the rows representing all the individual scan positions and the columns the energy channels. The disadvantage to such a sparse representation is that querrying/slicing/manipulation of the data is not straightforward at all, and without metadata such as the number of frames and scan dimensions the dataset is ambiguous. So spectrum stream data as handled by TEMMETA must always be accompanied by metadata.
Let's load a SpectrumStream
into memory. We again use the get_dataset
command to load in a SpectrumStream this time (see beginning of this file). You can verify that this object will be of the type SpectrumStream
.
You may get a warning that there is no FrameLocationTable
. Datasets with newer versions of Velox include a long array stating at which index in stream each frame begint. In older datasets this is not the case. If there is no FrameLocationTable, a limited number of functions does not work, and the number of frames must be infered from looking at other datasets in the file. However it has no further bearing on working with the dataset
specstr = emd1.get_dataset("SpectrumStream", "f5a4ba0965a5444b8c46cc420cf7fef0")
ERROR:temmeta.data_io:No FrameLocationTable found in f5a4ba0965a5444b8c46cc420cf7fef0
You can check the type of the data stored in specstr
specstr.data
<15728640x4096 sparse matrix of type '<class 'numpy.uint16'>' with 202136 stored elements in Compressed Sparse Row format>
You can querry an individual frame and this will be converted to SpectrumMap
object; this will store a full 3D matrix representation. Because only 1 frame has been selected, it will most likely be very empty and you will likely rarely use this command
frm123 = specstr.get_frame(123)
type(frm123)
temmeta.data_io.SpectrumMap
Alternatively, you can create new SpectrumStream
objects with select_frames
and exclude_frames
. If something happened during the acquisition that you want to select (or neglect) then you should consider these functions.
print(f"Number of frames in original: {specstr.frames}")
frams = [2, 4, 56]
spec_select = specstr.select_frames(frams)
print(f"Number of frames in spec_select: {spec_select.frames}")
spec_excl = specstr.exclude_frames(frams)
print(f"Number of frames in spec_excl: {spec_excl.frames}")
Number of frames in original: 240 Number of frames in spec_select: 3 Number of frames in spec_excl: 237
You can export the SpectrumStream
in two ways: exporting the individual frames as individual sparse matrix files in a folder, or as one sparse matrix file. The resulting files can not be read by other programs, they mainly serve to be shared with others. Along with the export, a metadata json file is stored by default. This file must be present in the same folder as the npz file, otherwise the npz file can not be re-imported properly.
Here we demonstrate exporting and re-importing all the frames as individual npz files
export_path = "./Test"
specstr.export_streamframes(export_path)
specstr_import = dio.import_files_to_spectrumstream(export_path)
Here we demonstrate exporting and re-importing the entire stream as one npz file
export_path2 = "./Test2/teststreamframe"
specstr.export_data(export_path2)
specstr_import2 = dio.import_file_to_spectrumstream("./Test2/teststreamframe.npz")
Once the appropriate frames have been selected and aligned in the stream, we can do 2 main operations: create a summed spectrum or create a spectrumstream. Here we look at the spectrum, as it's quite similar to an intensity profile. Essentially we add all rows in the sparse matrix which gives us our spectrum; this is an extremely efficient operation. We calculate the spectrum and return a Spectrum
object with the spectrum
attribute
spec = specstr.spectrum
We can use the plot
command for visulatization. Again, we can edit the properties of the plot by catching the objects returned by the plot command
ax, spec_plt = spec.plot()
ax.set_ylim(0,1500)
ax.set_xlim(0, 20)
spec_plt[0].set_label("Total spectrum")
spec_plt[0].set_color("red")
ax.legend()
<matplotlib.legend.Legend at 0x16f6a24c0>
An additional function to the Spectrum
object is that you can find peaks using the calc_peaks
method. You should provide some arguments in the pf_props
dictionary, which is passed directly to c. Please familiarize yourself with the options here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html . The function will return a table of peaks with columns channel index, energy, and peak height. In addition, it will return all the additional properties of the peaks, calculated by scip.signal.find_peaks
.
spec.calc_peaks(pf_props={"height": 30, "width": 10})
( Channel Energy (keV) Intensity (a.u.) 0 148.0 0.264302 751.0 1 202.0 0.534302 693.0 2 237.0 0.709302 876.0 3 480.0 1.924302 71.0 4 530.0 2.174302 1017.0 5 1373.0 6.389302 1227.0 6 1507.0 7.059302 170.0 7 1700.0 8.024302 46.0 8 3403.0 16.539302 172.0, {'peak_heights': array([ 751., 693., 876., 71., 1017., 1227., 170., 46., 172.]), 'prominences': array([ 675., 478., 859., 38., 1015., 1227., 168., 46., 172.]), 'left_bases': array([ 118, 174, 118, 438, 118, 118, 1438, 1586, 3338]), 'right_bases': array([ 174, 212, 340, 500, 973, 1586, 1586, 1840, 3471]), 'widths': array([12.94681677, 10.83 , 15.8115722 , 12.27777778, 25.84555556, 26.79081633, 28.07857143, 18.25 , 38.22380952]), 'width_heights': array([413.5, 454. , 446.5, 52. , 509.5, 613.5, 86. , 23. , 86. ]), 'left_ips': array([ 143.09782609, 194.63 , 228.83974359, 470.5 , 520.21 , 1361.56632653, 1491.85 , 1693.5 , 3386.64285714]), 'right_ips': array([ 156.04464286, 205.46 , 244.65131579, 482.77777778, 546.05555556, 1388.35714286, 1519.92857143, 1711.75 , 3424.86666667])})
The basic peaks table is stored in the peaks
attribute. It is a pandas Dataframe, which you can export to excel if necessary with to_excel
(see documentation).
spec.peaks
Channel | Energy (keV) | Intensity (a.u.) | |
---|---|---|---|
0 | 148.0 | 0.264302 | 751.0 |
1 | 202.0 | 0.534302 | 693.0 |
2 | 237.0 | 0.709302 | 876.0 |
3 | 480.0 | 1.924302 | 71.0 |
4 | 530.0 | 2.174302 | 1017.0 |
5 | 1373.0 | 6.389302 | 1227.0 |
6 | 1507.0 | 7.059302 | 170.0 |
7 | 1700.0 | 8.024302 | 46.0 |
8 | 3403.0 | 16.539302 | 172.0 |
Once you have calculated the peaks you can also plot them by setting the show_peaks
argument to True
in plot
. It will plot the peaks as points as well as a red interval with w
the number of channels in this interval
ax, spec_plt = spec.plot(show_peaks=True, w=30)
ax.set_ylim(0,1500)
ax.set_xlim(0, 20)
(0.0, 20.0)
We can export a spectrum similar to an intensity profile, to an excel file or to a hyperspy object. With hyperspy, many more tools are available for analyzing/manipulating EDX and EELS spectra. The conversion is just to a simple signal 1D objecta and metadata is lost. You have to manually convert it to an edx or eels dataset and set the important metadata (see http://hyperspy.org/hyperspy-doc/current/user_guide/eds.html). In future versions we can consider automatic recognition of EELS vs EDX and automatic conversion.
spec.to_excel("./spectrum.xlsx")
spec.to_hspy(filepath = "./spectrum")
<Signal1D, title: , dimensions: (|4096)>
You can create a summed spectrum map with the attribute spectrum_map
smap = specstr.spectrum_map
As with the spectrum stream you can also still access the spectrum from the spectrum map. Changes to the spectrum map will also be implemented in the spectrum
spec_2 = smap.spectrum
spec_2.plot()
(<AxesSubplot:xlabel='Energy (keV)', ylabel='Intensity (a.u.)'>, [<matplotlib.lines.Line2D at 0x15be56e80>])
You can visualize the spectrum map with plot_interactive
. This will show an image corresponding to the spectrumstream integrated over a certain energy window. The energy slider sets the middle of this window, the width determines the width of the window. Moving the energy slider allows you to explore different parts of the spectrum, to help you the spectrum is plot on the slider. I haven't found a way to make the slider taller unfortunately. I recommend clicking to a position in the energy slider and not actually sliding: this can block the widget after a while. The max slider controlls the color mapping.
smap.plot_interactive()
You can perform basic operations on the spectrum, such as rebinning and cropping. Rebinning the spatial dimension is a different command from rebinning the energy dimension. Note that rebinning works differently from rebinning images: here it works in an additive way! Neighboring pixels are simply summed together. Hence we only accept integer divisors of the dimension sizes as argument here. If you need more advanced binning with linear interpolation consider using hyperspy for these operations.
smap.rebin(2, inplace=True)
smap.crop(20, 0, 60, 70, inplace=True)
smap.rebin_energy(4, inplace=True)
smap.data.shape
(1024, 70, 40)
smap.plot_interactive()
You can extract certain "peak" images from the map with the integrate_to_image
command where you provide the energy and the window width in the correct units. Returned is a GeneralImage object that you can work with as before. Here we use the peaks calculated earlier to select the energy position
spec_img1 = smap.integrate_to_image(spec.peaks.iloc[5, 1], 0.4)
spec_img2 = smap.integrate_to_image(spec.peaks.iloc[4, 1], 0.4)
spec_img1.plot(imshow_kwargs={"cmap": "jet"})
(<matplotlib.axes._axes.Axes at 0x14a509dc0>, <matplotlib.image.AxesImage at 0x14aafe910>)
You can convert the spectrum map also to a hyperspy dataset and continue working with it there
It is also possible to extract line spectra from spectrum maps, analogous to extracting intensity profiles in images. A wedge is taken out of the data cube along the user defined direction. So the data in the line spectrum is 2-dimensional: length x channels.
The way it is currently implemented, calculation of such a profile can be quite slow with large data cubes (a few minutes even for $256 \times 256 \times 4096$). In any case, you will get very few counts in a line extracted from a map. Therefore, I highly recommend you do a binning first on the map in both the spatial and energy directions
smap = specstr.spectrum_map
smap.rebin(2, inplace=True)
smap.rebin_energy(8, inplace=True)
smap.data.shape
(512, 128, 128)
To demonstrate where we will make the line profile, we first plot an image from one of the peaks we saw before in the spectrum. I just selected two random coordinates for the green arrow and will only integrate over a width of 5 pixels. For the magenta arrow we take one directly perpendicular to the linear features and will integrate over a wide range
im_bigpeak = smap.integrate_to_image(6.51, 0.5)
ax, impl = im_bigpeak.plot()
impl.set_clim(0, 8)
impl.set_cmap("hot")
x1, y1 = (35, 103)
x2, y2 = (70, 13)
ax.arrow(x1, y1, x2-x1, y2-y1, color = (0., 1., 0.), head_width=5)
x1b, y1b = (64, 113)
x2b, y2b = (64, 6)
ax.arrow(x1b, y1b, x2b-x1b, y2b-y1b, color = (1., 0., 1.), head_width=5)
<matplotlib.patches.FancyArrow at 0x14c2037f0>
The line profile is calculated with the line_spectrum
command. It works in the same way as getting a line profile
reload(dio)
<module 'temmeta.data_io' from '/Users/nielscautaerts/Documents/PythonProjects/TEMMETA/temmeta/data_io.py'>
slp = smap.line_spectrum(x1, y1, x2, y2, w=3)
slp2 = smap.line_spectrum(x1b, y1b, x2b, y2b, w=50)
DEBUG:temmeta.data_io:Shape of line profile data: (512, 92) DEBUG:temmeta.data_io:Shape of line profile data: (512, 107)
We can visualize the spectrum with plot_interactive
. With the sliders you can select an energy and an energy integration width. With slp2
, we can clearly see the evidence of the atomic rows. Select the peak around 2.23 KeV and you will see maxima where the dark bands are. Selecting 6.51 will of course show minima there.
slp2.plot_interactive()
To create proper profile objects for plotting at a particular energy, you use the get_profile
command and provide the energy and integration window. This will return a Profile object which is identical to intensity profiles in images. You can then make plots as was shown early on in this document.
prof_651 = slp2.get_profile(6.51, 1)
prof_651.plot()
(<AxesSubplot:xlabel='x (m)', ylabel='Intensity (a.u.)'>, [<matplotlib.lines.Line2D at 0x15b64b280>])
prof_651.history
* Integrated channels 162-187 (6.51 p.m. 1 keV) | * Spectrum line profile from (64, 113) to (64, 6). Integration width: 50 | * Rebinned channel axis by factor 8 | * Rebinned scan_x and scan_y axis by factor 2 | * Sum of all frames | * Extracted dataset SpectrumStream/f5a4ba0965a5444b8c46cc420cf7fef0 from EMD file stem-edx.emd
Finaly, you can export a line profile to a hyperspy dataset as with most TEMMETA objects
slp_hs = slp.to_hspy()
I hope TEMMETA will be useful for analyzing your TEM data.
- Niels Cautaerts